A weighted message-passing algorithm to estimate volume-related properties of 

random polytopes 



Franccsc Font-Clos 

Centre de Recerca Matematica, Edifici C, Campus UAB, E-08193 Bellaterra (Barcelona), Spain 

Francesco Alessandro Massucci and Isaac Perez Castillo 
Department of Mathematics, King's College London, WC2R 2LS United Kingdom 

In this letter, we introduce a novel message-passing algorithm for a class of problems which can 
be mathematically understood as estimating volume-related properties of random polytopes. Unlike 
the usual approach consisting in approximating the real-valued cavity marginal distributions by a 
few parameters, we propose a weighted message-passing algorithm to deal with the entire function. 
Various alternatives of how to implement our approach are discussed and numerical results for 
random polytopes are compared with results using the Hit-and-Run algorithm. 



As it is well-known, old and new interesting problems 
such as Gardner's optimal capacity problem for contin- 
uous synaptic couplings von Neumann's problem of 
linear economics estimation of reaction fluxes from 
mass-balance equations l3| , or reconstruction techniques 
in compressed sensing [J|, can be related, in one way or 
another, to estimating volume-related properties of ran- 
dom polytopes. In the so-called H-representation, a poly- 
tope is defined as the set of points x = {xi, . . . ,xn) e 
D CR^ encapsulated by M hyperplanes {(|^, 7'')}^ii, 
with = (^f,...,^^) the normal vector of the ^- 
hyperplane. From all possible questions related to poly- 
topes in H-representation, we are particularly interested 
in that of the volume V and its projection onto each axis 
relative to V or, in order words, the marginal pdf Pi{xi). 
These two quantities can mathematically be written as: 
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where G(a;) is the Heaviside step function and x\i denotes 
the vector x without component i. 

Finding efficient methods to obtain reliable estimates 
for the marginals Pi{xi) has been, and still is, the main 
objective of past and present research. In the context of 
metabolic networks, the area we will focus on to present 
our work, these methods can roughly be divided into ei- 
ther theoretical approaches of some sort, or Monte Carlo 
simulations. In the first class of methods, we find Flux 
Balance Analysis (FBA) Q, which consists in approxi- 
mating the whole volume V of plausible solutions of the 
mass-balance equations by a single point, which is se- 
lected by optimising an objective function [l6[. FBA pro- 
vides reliable results when applied to elementary micro- 
organisms in simple milieu conditions (see, for instance, 
0, IE]), but the collapsing of to a single point is too 
drastic for more complex micro-organisms. To explore 



the whole volume V one usually relies on Monte Carlo 
simulations, as the one introduced in Q, or the Hit-and- 
Run algorithm 0. These two algorithms provide a uni- 
form s amp ling of V, but they are restricted to small net- 
works |17| ■ For larger networks the MinOver+ algorithm 
has shown to be very promising , with the sampling 

becoming more uniform the larger the network is. 

Here we want to re-explore the possibility of using 
message-passing algorithms to estimate Pi{xi). Message- 
passing equations, that is cavity equations, were already 
derived in Q for the von Neumann problem, and in 11 1 
for mass-balance problem. As the dynamical variables 
are continuous, it was thought that having a message- 
passing algorithm for the exact marginals was a daunt- 
ing task |8| , even though fast Fourier transform has been 
nicely applied in some cases |11|. 

Let us see how it is possible to introduce a weighted 
message-passing algorithm for the entire cavity density 
marginals. Having applications in the area of metabolism 
in mind, we explain our methodology in the context 
of the von Neumann model applied to metabolic net- 
works SElii- The application to similar problems is 
straightforward. In here, a metabolic network is a sys- 
tem of i ~ 1,...,N coupled chemical reactions which 
produces and consumes fi = 1,...,M metabolites and 
exchanges them with the environment. Following von 
Neumann's work [l, as some of the the input will be 
used to produce output, the ratio of input produced to 
output consumed cannot be larger than a global growth 
rate p, that is 
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where A — (af ) and B ~ (b^) are the input and output 
matrices (e.g. stoichiometric coefficients), Xi is the rate 
of chemical reaction i, and 7 = (7^, . . . , 7*''^) is the vec- 
tor of exchanged rates (associated to e.g. transport of 
metabolites). Given a metabolic network, that is, given 
a set of stoichiometric coefficients and exchanged rates. 
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we want to find the volume of solutions V to the set of 
inequalities (|3]) at the optimal growth rate. As we do not 
want to get distracted by unnecessary complicacies, we 
denote = (p) = — paf , ignoring the dependence 
on p and calling ^ = (^l") the stoichiometric matrix. 
In practical situations, the stoichiometric matrix is gen- 
erally diluted and it is not unreasonable to consider the 
problem of the volume for diluted random polytopes, 
meaning that the corresponding bipartite graph associ- 
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ated to the stoichiometric matrix £ is tree-like for both 
i-nodes and ^-nodes. The volume of solutions V is given 
by ([1]) . with D being the domain region of the reaction 
rates. This domain is usually determined by the biochem- 
istry of each reaction, that is, reaction rates are usually 
bounded Xi G [x™"^, xf"^"]. For sake of simplicity we drop 
the domain from the notation. Using the cavity method, 
we obtain the following set of coupled cavity equations 
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where mj^ and are normalisation constants and the 
notation e.g. j E 9/i means considering all reactions 
which participate in the production or consumption of 
metabolite fx. Here s = (xi, . . . , xjv, zi, . . . , zm), where 
the z— variables have been introduced to write the step 
functions as integrated Dirac deltas, while the matrix 
rf corresponds to the augmented stoichiometric matrix 
resulting from this transformation. Although this is not 
a necessary mathematical step, it will be useful for the 
subsequent discussion. 
The actual marginals Pi{xi) are given by 

^'(^'■) = ^ n "'-^(^')' * = l:---,^. (6) 

while the meaning of m^*''(xi) in eq. ^ is the follow- 
ing: assuming that the marginals {P^^^^ {x^)} ifz^^jy^i when 



reaction i has been removed are known, the probability 
that, by fixing the value at node i to Xi the inequality /i 
will be fulfilled, is precisely m^*''(xj;). 
At this point, one would be tempted to use the Dirac 
delta in eq. Q in an iteration process with assign- 
ments (7'' - J2jedt^\i'njSj)/T]'^ Si as it is usually 
done in the standard disordered-averaged cavity equa- 
tions [3]. There are two differences, though: (i) the 
set of equations Q is still on an instance; (ii) due 
to the integration domain, not all possible assignments 
(7^ ^ J2jedtJ.\i''lj ^ ^^'^ allowed. Fortunately, 

the acceptance region in eq. ^ for m\^\si) is easy to 
determine as it corresponds to restricting the domain D 
by the hyperplane /i. Within the acceptance region, we 
can rewrite equation Q as follows: 
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Here fc^'^ ~ \dfi\i\ and £j e 9/i \ t for j = 1, . . . , fc^*^ The new expression ([7]) invites us to propose the following 
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updating method to solve numerically this equation: 

1. Draw with probability P^'^\s£-^\^'^), draw se^ 

. . ., draw Si ^.j 



with probability Pff {si^\si^;^'^) 
with probability p/^' (si |s£i, . 



2. Assign ^ (7^ - J2j^a^\,VjSj)/Vr with weight 
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With this weighted iteration in mind, it is possible to 
solve numerically the set of equations Q and (O directly. 
We propose two methods: (i) method of histograms. 
Here one represents {Pi'^\sij)} with histograms. Firstly, 
we use the weighted iteration method in equation ^ to 
obtain a population of size TV of pairs {{si,a, ^i^a )}c[=i 
from which to construct a histogram for mlf' [si) = 
(1/ Ell Ell - Then use equa- 

tion ([5]) to recalculate the histograms of {F^^^^(s(,J} from 

the histograms of {m^*^(si)}; (ii) method of the weighted 
populations. It is possible to circumvent the step of using 
and constructing histograms during the iteration process. 
One simply works directly with a set of equations for the 
cavity marginals {Pi^\st,.)}, by combining eqs. and 
([5]) • Each of these functions is represented by a weighted 
population {(sj,Q, rj^2)}li of size TV. Note that in this 
case, one Dirac delta appearing in the combined expres- 
sion gives the assignment for Sj while at the same time, 
locks the other Dirac deltas, providing an extra weight. 

To assess the validity of our approach we have run 
some numerical checks and compared them with Monte 
Carlo simulations using the Hit-and-Run algorithm. We 
entirely focus our numerics in the methods of histograms 
and postpone a more detailed analysis on the method 
of weighted populations for another occasion (T^]. We 
have created a random metabolic network of = 25 
reactions and M = 10 metabolites. The in- and out- 
degree distributions for reaction and metabolite nodes 
are chosen to be Poisson with average degree 1.5 and 3.75, 
respectively. The stoichiometric coefficients as well as 
the exchange fluxes 7^ are chosen randomly, but ensuring 
that the volume of solutions is not empty. Finally, we 
assume that the domain D = [0, 1]^^. 

Results are summarised in fig. [1] which shows all 
marginals {Pi{xi)}'^li obtained by the weighted message- 
passing algorithm, together with the Monte Carlo simula- 
tions using the Hit-and-Run algorithm. The comparison 
is fairly excellent, considering that the cavity equations 
are based on the Bethe approximation, which does not 
generally hold for a random graph of finite size, particu- 
larly this small. 

Looking at the results in fig. [T] we also note that, 
although originally the domain for each reaction rate 
is the interval [0,1], some reaction rates seem to have 
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Figure 1: Result for the marginals Pi{xi) for i = 1,...,25 
(sliown in red and ordered from left to right and up to bottom) 
using the weighted iteration algorithm within the histogram 
method for a random metabolic network of A'' = 25 reactions 
and M = 10 metabolites. Our numerical results are com- 
pared with Monte Carlo simulations using the Hit-and-Run 
algorithm (green line). Finally, the solid blue circles are the 
results of the belief-propagation algorithm for the boundaries 
of the marginals {Pi[xi)~\. 



a smaller domain region. This can be due because ei- 
ther the polytope further constrains the original domain 
or rather because there exists a positive but very small 
probability. It is important to distinguish between what 
cannot happen and what is unlikely to happen, particu- 
larly in metabolism, where a perturbed network can be 
brought into a region of originally unlikely events (i.e. 
a disease). To shed some light into this matter, we no- 
tice that the set of equations (j4|) and ([5]) suggests a way 
to write a belief-propagation algorithm for the domains 
[flj, bi] 9 Xi of the marginals Pi{xi). (see [l3| for details). 
The results are also reported in fig. [T]as solid blue circles 
and indicate that, for this particular instance, the prob- 
ability in some regions for some variables is very small. 
To explore further into this matter and to access these 
regions, we have performed a very simple perturbation 
of the original metabolic network consisting in changing 
the domain D by [0, bY^ for b = 1.0, 0.60, 0.30, 0.12, and 
h = 0.10. Results are shown in fig. [2] for variable 18, 
where we can see that not only the probability increases, 
but also that the shape of the pdf changes as well in a 
non-trivial way. 

It is interesting to keep in mind that since a marginal 
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Figure 2: Results for the marginal Pi{xi) for i = 18 and 
by varying the hypercube [0,fe]^^ for h = 1.0,0.60,0.30,0.12, 
and b = 0.10. For a visually easier comparison, the hor- 
izontal and vertical axes have been rescaled as Xi/b and 
fi{x) = Pi{xi)/Mi with Mi = max^ (z[o^t] Pi (x), respectively. 



to deal with the exact cavity marginals, rather than ap- 
proximating these with a set of parameters. We have 
introduced two ways of implementing the method and 
discussed here one of them, the method of histograms. 
This method has admittedly the drawback of having to 
approximate the functions with histograms, but avoids 
the rejection region characteristic of the problem. 

Among the various research lines were are currently 
considering, we would like to mention two. Firstly, we 
are currently exploring the implementation of the algo- 
rithm by the method of weighted populations [l3l |. Like 
in the standard replica symmetric equations in the en- 
semble, the cavity marginals are naturally represented by 
a population of random variables. In this way, we avoid 
having to introduce any discretization of the marginals. 
Secondly, we are extending this approach to the more in- 
teresting case of the double von Neumann problem fl5| . 
which is able to incorporate information on the energy 
balance of the chemical reactions. 



Pi{xi) corresponds to projecting the volume of the poly- 
tope onto the x^-axis, restrictions in the polytope may 
affect both boundaries [a, 6] 9 x^. To follow track of 
these changes the message-passing equations for the do- 
main [l^l are most convenient. As an example, in fig. 
Owe have plotted how the domain [a, b] of the marginal 
Pi{xi) for variable 22 changes as a in the original domain 
[a, 1]^'^ of the polytope is varied. 



T 1 1 r 




0.25 



a 



Figure 3: Behaviour of domains [a, 6] 9 Xi of marginal Pi{xi) 
variable i = 22 as the end point a is varied for a perturbed 
domain [a, 1]^^. The inset captures a more detailed behaviour 
of the region where the end points meet. 

In this letter we have presented a weighted message- 
passing algorithm to deal with a class of problems which 
consists in counting solutions of a set of equations (ei- 
ther equalities or inequalities). The method is designed 
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